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Abstract 

We performed the detailed NLO analysis of the combined CCFR xF% and F2 struc- 
ture functions data and extracted the value of a s , parameters of distributions and 
higher-twist (HT) terms using the direct solution of the DGLAP equation. The value 
of a s {M z ) = 0.1222 ± 0.0048(exp) ± 0.0040 (t/ieor) was obtained. The result has larger 
central value and errors, than the original result of the CCFR collaboration, in view of 
the incorporation into the fits of the HT terms as the free model independent param- 
eters. The x-shapes of the HT contributions to xF% and F2 are in agreement with the 
results of other model-independent extractions and are in qualitative agreement with 
the predictions of the infrared renormalon model. We also argue that the low x CCFR 
data might have the defects, since their inclusion into the fits led to the following low 
x-behaviour of the gluon distribution xG(x,9 GeV 2 ) ~ x o.092±o.oo73^ - n CO ntradiction 
with the results of its extraction from low x HERA data. 
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1. The study of the possibility to separate power suppressed terms (namely higher- 
twist (HT) effects) from the perturbation theory logarithmic corrections in the analysis of 
scaling violation of the deep-inelastic scattering (DIS) processes has a rather long history 
(see e.g. Refs. |l], 0, |3| and Ref.|| for the review). In the recent years the interest to this 
problem was renewed, mainly due to the consideration of the possibility to model HT terms in 
different processes using the infrared- renormalon (IRR) technique (see e.g. Refs.[||]-[10|,[11] 
and especially Ref.0| for the review). 

On the other hand the experimentalists are improving the precision of their data and 
are achieving, sometimes, percent level of accuracy. For example the data on xF% and F 2 
from the most precise vN DIS experiment, performed at Tevatron by CCFR collaboration, 
recently appeared [|l|, [lij]. The CCFR data on xF 3 were analysed in Ref.fL5| in the LO 
and and with inclusion of the NLO and approximate next-to-next-to-leading order (NNLO) 
corrections. For the latest, the NNLO QCD corrections to the coefficient function |16[ were 
taken into account. The NNLO corrections to the anomalous dimensions of a limited set 
of even non-singlet moments ]T7|] were also taken into account. The NNLO corrections to 
the anomalous dimensions of odd moments, which are not still explicitly calculated, were 
obtained using smooth interpolation procedure proposed in Ref.fl8|| and improved in Ref . |19 1 . 



The aim of the work of Ref.|l5|] was to make an attempt of the first NNLO determination of 
a s (Mz) from DIS and to extract the HT terms from the data on xF 3 within the framework 
of the IRR-model and also by the model-independent way, similarly to the analysis of 
the combined SLAC/BCDMS data ||20|| , performed in the NLO approximation. Theoretical 
uncertainties of this analysis were further estimated in Refs. [pT], 22] at the N 3 LO using the 
method of Pade approximants. It was found in Refs. 



15, 21, 22 



that the inclusion of the 
NNLO corrections leads to the decrease of the HT contribution value, so that at the NNLO 
its x-shape lies closer to zero. 

In these analyses only statistical errors of data were taken into account. However, the 



systematic errors of the CCFR experiment are not small |LJ] and may even dominate in the 
determination of some parameters. In this paper we filled in this gap and performed the 
NLO analysis of the CCFR data with the help of QCD DGLAP evolution code, developed in 
Ref. P5fl . (Remind, that the analyses of Refs.[|T5|, |2T], [22|] were performed with the help of the 
Jacobi polynomial variant |^4[ of the DGLAP equation [p7| ). In addition, we included 

in our analysis the CCFR data on the singlet structure function F 2 . It should be stressed 
that the code p3| was tested using the procedure proposed in Ref. J28| and demonstrated 
the accuracy at the level of 0(0.1%) in the kinematic region covered by the analysed data. 
It was already applied for the nonsinglet DGLAP analysis of the combined SLAC/BCDMS 
SLAC data on F 2 f2| . 

2. Our fits were made in the NLO approximation within the MS factorization and 
renormalization schemes. The Q 2 dependence of the strong coupling constant a s was defined 
from the following equation 
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where (3q = 11 — (2/3)n/ and /3i = 102 — (38/3)n/. Note, that the explicit solution of Eq.(l) 
can be expressed through the Lambert function [fjQfl . However, we did not use in our work 
this explicit representation and solved Eq.(l) numerically. The effective number of flavours 
rif was chosen to be rif = 4 for Q 2 less than the definite scale Mf and increased to rif = 5 
at larger values of Q 2 keeping the continuity of a s [31]. The value of the effective matching 
scale M 5 was varied from M 5 = mj to M 5 = 6.5m;,. The last choice was advocated in Ref.|32| 
on the basis of the DIS sum rules consideration. The dependence of the results of the fits on 
the choice of the matching point gives one of the sources of theoretical uncertainties inherent 
to our analysis. 

The leading twist term xF^ T (x,Q) was obtained by direct integration of the DGLAP 



equation f27 
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Ots{Q) 
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where P qq (x) denotes the NLO splitting function, taken from Ref.[]34|. The function xF 3 is 
determined by the subsequent convolution with the NLO coefficient function C^ q (x) : 



xF^ T (x, Q) = f 1 dzQ 

JX 



X 

A/ 



NS, 



x/z,Q). 



(4) 



The boundary condition at the reference scale Q"l 



5 GeV 2 was chosen in the form analogous 



to the ones, used in Refs. [1141 ITS 
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and t]ns is the measure of the deviation of the Gross-Llewellyn Smith integral [35] from its 
quark-parton value 3. The expression for the xF 3 , which includes the HT contribution, looks 
as follows: 

H 3 (x) 



xF^ 1 (x,Q) = xFr 



LTTMC, 



x,Q) + 



Q 2 



(7) 



where F^ T,TMC (x, Q) is F^ T (x, Q) with the target mass correction [^6| applied. 



At the first stage of this work, in order to perform the cross-checks of the code against 
the results of Refs.[[L4], [13], we fitted the data on xF 3 in the kinematical region Q 2 > 5 GeV 2 , 
W 2 > 10 GeV 2 , x < 0.7 (the number of data points is NDP=86). We made three fits with 
various ways of taking into account the HT effects. The first fit with no HT, i.e. H 3 (x) = 0, 



was done to compare our results with the ones of Table 1 of Ref.|L5|], obtained using different 
method EI fE 



and different computer code. The second fit with the HT chosen as 
one-half of the IRR model predictions [0, i.e. 



H 3 (x) = A> 2 f 1 dzCi RR (z)-F 3 LT (x/z,Q), 

Jx Z 



where 



and At, 



(jIRR 



1 



+ 2(2 + z + 2z 2 ) - 55(1 - z) - 5'(1 - z) 



(8) 
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-0.1 GeV 2 , advocated for the first time in Ref. IHI- The aim of this fit was to 



compare its outcomes with the results of Table 7.9 of Ref.[]T4|, where the computer code 
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Table 1: The results of the fits to data on xF% with statistical errors only I) without HT-terms, 
II) with HT accounted as one-half of the IRR model predictions, III) with model independent 
HT-contributions. HP' m4) ' m8) are t h e va i ues Q f # 3 ( x ) at x = . ; q.2, 0.4, 0.6, 0.8. 





I 


II 


III 


X 2 /NDP 


88.5/86 


81.9/86 


70.3/86 


b 


0.789 ± 0.024 


0.786 ± 0.024 


0.805 ± 0.067 


c 


a no i nil 
4.UZ zt U.ll 


A fin 1 n 11 
4.UU ± U. 11 


A OA -J- n o~\ 
4.Z4 ± U.Z1 


/ 


29 ± 30 


26 ± 30 


61 ± 71 


Vns 


0.927 ±0.014 


0.949 ±0.014 


0.927 ±0.030 




0.1193 ±0.0025 


0.1219 ±0.0024 


0.1216 ±0.0066 


HP 






0.18 ±0.19 


hP 






-0.26 ±0.12 








-0.21 ±0.31 


HP 






0.11 ±0.26 


HP 






0.90 ±0.47 



written by Duke and Owens was used. In the third fit we used the model independent 
HT-expression, i.e. H 3 (x) parametrized at x = 0.,0.2, 0.4,0.6,0.8 with linear interpolation 
between these points. It was performed to compare our results with the ones presented in 
Table3ofRef. JTJ]. 

All results of these our fits are presented in Table 1. We observed a good agreement of our 
results on a s with the both referenced papers. However, in the case of the values of x-shape 
parameters we found the certain discrepancy with the results of Ref . ||1 5|| . For example, the 
value of 7, as presented in column I of Table 1, is 7 = 0.26 ± 0.30, meanwhile the analogous 
parameter in Ref.|l5| is 7 = 1.96 ± 0.36. At the same time our x-shape parameters are in 
agreement with the ones, extracted in Ref. [14] within errors. In addition, we made the fit, 
releasing parameter A 2 and obtained the value of A 2 = —0.12 ± 0.05 in agreement with the 
results of Ref. 
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3. The next step of our analysis was to take into account the point-to-point correlations 
of the data due to systematic errors which, as we mentioned above, can be crucial for the 
estimation of full experimental errors of the parameters (see in particular |2"9" |, where the 
value a s {Mz) = 0.1180 ± 0.0017 (stat + syst) was obtained as the result of the combined 
fits of the SLAC/BCDMS data with HT included). The systematic errors were taken into 
account analogously to the earlier works |23|, ^9|. The total number of the independent 
systematic errors sources for the analysed data is 18 and all of them were convoluted into 
the general correlation matrix, which was used for the construction of the minimized x 2 - The 
results of the fits to xF 3 data with the model independent HT and with the systematic errors 
taken into account are presented in the first column of Table 2. One can see that the account 
of systematic errors leads to the significant increase of the experimental uncertainties of the 
extracted HT contributions (compare the first column of Table 2 with the third column of 
Table 1). In addition, the central values of the HT parameters moved. However, even in 
this case there is definite agreement with the results of HT-behaviour of Ref. , obtained 
in NLO. Moreover, these results do not contradict to the IRR-model prediction of Ref. 0, 
since releasing A 2 we obtained A 2 = —0.10 ± 0.09. 

Trying to minimize the errors of the parameters we added to the analysis the CCFR data 
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Table 2: The results of the fits with the account of systematic errors and model independent 
HT-effects. #(°,)>( 2 )>( 4 )>( 6 )>( 8 ) are the values of H ^ and H ^ x = 0., 0.2, 0.4, 0.6, 0.8. I) 

xF 3 with the cut Q 2 > 5 GeV 2 , Q 2 = 5 GeV 2 II) xF 3 kF 2 with the cut Q 2 > 5 GeV 2 , 
Ql = g GeV 2 III) xF 3 &F 2 with the cut Q 2 > 1 GeV 2 , Q 2 = 9 GeV 2 . 





T 
1 


TT 
11 


TTT 
in 


y 2 /NDP 


55.7/86 


154.9/172 


204.2/220 


^NS 


0.797 ±0.076 


0.800 ±0.016 


0.782 ±0.014 


CNS 


4.24 ±0.21 


4.060 ±0.068 


4.131 ±0.056 


7 


0.75 ±0.79 


0. 


0. 




0.945 ± 0.043 


0.922 ±0.027 


0.920 ±0.025 


a s (M z ) 


0.1269 ±0.0065 


0.1248 ± 0.0048 


0.1131 ±0.0045 






0.1785 ±0.0077 


0.1796 ±0.0065 






0. 


-0.034 ±0.023 


C5 




8.37 ±0.21 


8.00 ±0.29 






0. 


0.092 ±0.073 


CG 




7.5 ±2.6 


11.50 ±0.90 


?7g 




0.69 ±0.35 


1.08 ±0.19 


flf 




-0.23 ±0.56 


0.09 ±0.11 


lf< 2 > 




-0.28 ±0.18 


-0.239 ±0.094 


1#> 




-0.14 ±0.18 


0.17 ±0.13 


tff 




-0.03 ±0.13 


0.204 ±0.097 


iff 




0.21 ±0.18 


0.14 ±0.18 




0.28 ±0.21 


0.34 ±0.11 


0.115 ±0.031 




-0.22 ±0.19 


-0.24 ±0.16 


-0.16 ±0.16 




-0.42 ±0.35 


-0.22 ±0.22 


0.28 ±0.19 




-0.09 ±0.28 


-0.05 ±0.17 


0.19 ±0.15 




1.21 ±0.50 


0.89 ± 0.44 


0.88 ±0.44 
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Figure 1: The high-twist contribution to F 2 . 



for the structure function F 2 . To perform the QCD evolution of F 2 one is to involve into the 
analysis the singlet and gluon distributions: 



iLT , 



x,Q) 



dz 



The distributions q F (x, Q) and G(x, Q) were obtained by integrating the system 



x 



C 2 , q (z)-(q NS (x/z, Q) + q PS (x/z, Q)) + C 2 , G (z)-G(x/ z, Q) 



x 



dxq p 


_ a,{Q) 


d\nQ 


IT 


dxG 


a s (Q) 


d\nQ 


7T 



dz 



dz 



x 



X 



P q y q b (z)-q yb (x/z,Q) + P qG (z)-G(x/z : Q) 
HH z z 



X 



X 



P Gq (z)-q^(x/z, Q) + P GG (z)-G(x/z, Q) 



with the boundary conditions 



xq 



PS 



(x,Qo)=Vsx bs (l-x) cs /A s , 



xG{x,Q ) =r] G x bG {l-x) CG /A c 



where 



A G 



x bs (l — x)° s dx, 

1- < xQ(x) > 
Jo x b( 3(l - x) c Gdx' 



(10) 



(11) 
(12) 

(13) 
(14) 

(15) 
(16) 



and < xQ(x) > is the total momentum carried by quarks. 

In order to provide the straightforward way for the comparison of our results with the 
analysis of Ref . p3| , the initial reference scale Qq=9 GeV 2 was chosen. In addition to the 
point-to-point correlation of the data due to systematic errors, the statistical correlations 
between F 2 and xF 3 were also taken into account. Performing the trial fits we convinced 
that adding the factor (1 + jx) to the reference expressions for the the gluon and singlet 
distributions do not improve the quality of the fit. Also we fixed parameters 7tvs, bs and b G 
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Figure 2: The high-twist contribution to xF% 



at zero because this increased the value of y 2 by few units only while x 2 /NDP remained 
less than unity. The HT contribution to F 2 was accounted analogously to xF 3 as: 



iHT, 



x,Q) 



F 2 (x, Q) + 



where H2(x) was parametrized in the model independent form. The results of the fits of the 
parameters H 2 (x) and H%{x) are presented in the second column of Table 2 and are depicted 
in Fig.l and Fig.2. One can note, that, comparing with the fit to xF 3 data only, the HT 
parameters errors are decreasing. Within the errors, the parameters, which describe the 



boundary distributions, are compatible with the outcomes of the similar fits from Ref.jn . 
The coefficients of H^{x) are in agreement with the NLO results of Ref. fl5| and the behaviour 
of H 2 (x) qualitatively reproduce the HT contribution to F 2 , obtained from the combined fits 
of the SLAC/BCDMS data on F 2 , performed in Refs. 



When the matching scale M5 was changed from m& to 6.5mb, the value of a s (Mz) shifted 
down by 0.0052 and then we ascribe to a s (Mz) the theoretical error of 0.0026 due to uncer- 
tainty of b-quark threshold matching. This uncertainty is in agreement with the results of 
the NLO Jacobi-polynomial fits of the CCFR data obtained within so-called spline MS pre- 
scription ]37| . One more source of the theoretical uncertainty, which is due to the truncation 
of higher QCD orders, was evaluated following the way, proposed in Ref. [p0| . In accor- 
dance with their procedure one can introduce renormalization scale k R into QCD evolution 
equations in the way, illustrated on the example of NS evolution: 



dxq 



NS 



dlnQ 



71 Jx I yy 



+ 



Ois{k R Q) r 



2tt 



z) + 



P g N q S ' W (z) + (3 P^\z) \n{k R ) \-q NS (x/z, Q) 



(17) 



where P NS ^°) and p NS '^-) denote the LO and the NLO parts of the splitting function P NS . 
The dependence of the results on k R would signal an incomplete account of the perturbation 
theory effects. The shift of a s (Mz) resulting from the reasonable variation of k R leads 
to an additional error of over 0.003 due to the renormalization scale uncertainty. Taking 
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Qq = 20 GeV 2 as the initial scale, we have checked that our results are rather stable to the 
variation of the factorization point. 

The NLO value we fare presenting as the main result is thus 

a s (M z ) = 0.1222 ± 0.0048(stat + syst) ± 0.0040(thresh + ren.scale) (18) 

It differs a bit from the NLO value a s {M z ) = 0.119 ± 0.002(stat + syst) ± 0. 00 A(theory), 



obtained in the CCFR analysis fL3| . The increase of the experimental errors is due to the 
fact that while CCFR group used model-dependent form of the HT contributions, we are 
considering them as the additional free parameters and are extracting them from the fits. 

In order to try to decrease further the errors, we repeated the fits of the combined xF% 
and F 2 data, using the less stringent cut Q 2 > 1 GeV 2 . The obtained results are presented in 
the third column of Table 2. In these fits the parameters b s and be were released since their 
values turned out to be statistically different from zero. We found, that the values of a s and 
be are correlated (the correlation coefficient is equal to -0.65). When we fixed be = 0, the 
value a s (M z ) = 0.1172 ± 0.0029 was obtained, and when we kept be as the free parameter, 
we obtained low value of a s [Mz) — 0.1131 ±0.0045. The analogous effect of correlations was 
observed for the fit with the cut Q 2 > 5 GeV 2 , although with less statistical significance. 
It should be underlined, that when we released bo in the fit with the cut Q 2 > 1 GeV 2 we 
faced another problem: its value turned out to be 

b G = 0.092 ±0.073, (19) 

which is in the evident contradiction with the results, obtained in the analysis of HERA data 
(for example the combined analysis of DIS from HERA and CERN-SPS data results in the 
value be = —0.267 ± 0.043 f2"3| , while in the framework of MRST parametrization the value 



be = — 1.08|38| was obtained). This problem might be related to the well-known discrepancy 
between CCFR and NMC/BCDMS data at small x. 
Conclusion 

In conclusion, we would like to stress that in order to perform similar analysis at the 
NNLO level it is necessary to calculate the yet unknown Altarelli-Parisi kernels to the corre- 
sponding DGLAP equations. Therefore, we are unable to obtain the results, similar to the 



NNLO ones of Refs.f^, [22|| . We hope that future progress of theoretical calculations will 
allow us to generalize our results to the NNLO approximation. 
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